#include "StVenKir.hpp"

// variational formulation language
//#include <feel/feelvf/vf.hpp>

#include <svk_cl.cpp>

namespace Feel
{


    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    StVenantKirchhoff<Dim, Order, Entity>::updateJacobian( const vector_ptrtype& X,
                                                           sparse_matrix_ptrtype& J,
                                                           vector_ptrtype& R)
    {
#if 0
        using namespace Feel::vf;

        if (M_LookComment) std::cout << "[StVenantKirchhoff] : updateJacobian start\n";

        this->changeRepository();

        //--------------------------------------------------------------------------------------------------//

        mesh_ptrtype mesh = M_Xh->mesh();

        element_struct_type u( M_Xh, "u" ); u = *X;
        element_struct_type v( M_Xh, "v" );

        //--------------------------------------------------------------------------------------------------//

        //coeff ...
        //   double lambda = this->vm()["rho"].template as<double>();
        //  double nu = this->vm()["rho"].template as<double>();


        //Matrice Identitee
        auto Id = oneX()*trans(oneX()) + oneY()*trans(oneY());

        auto Fv = Id + gradv(u);
        auto F = Id + grad(u);
        auto dF = gradt(u);
        auto dFF= grad(u);

        auto Ev = 0.5*(gradv(u)+trans(gradv(u)) ) + 0.5*trans(gradv(u))*gradv(u);
        auto E = 0.5*(grad(u)+trans(grad(u)) ) + 0.5*trans(grad(u))*grad(u);
        auto dE = 0.5*(trans(Fv)*dF + trans(dF)*Fv);
        auto dEE = 0.5*(trans(Fv)*dFF + trans(dFF)*Fv);

        auto Sv = lambda*trace(Ev)*Id + 2*mu*Ev;
        auto S = lambda*trace(E)*Id + 2*mu*E;
        auto dS = lambda*trace(dE)*Id + 2*mu*dE;
        auto dSS = lambda*trace(dEE)*Id + 2*mu*dEE;

        //--------------------------------------------------------------------------------------------------//
        form2( M_Xh, M_Xh, J, _init=true )  =
            integrate (elements(mesh), _Q<4*Order>(),
                       trace( (dF*Sv + Fv*dS)*trans(grad(v)) ) );

        //Attention les termes de symetrisation ne sont pas necessaire mais accelere la convergence
        // Prendre le parametre de stabilisation d'environ 1e6!!!
        double penalisation_bc=this->application()->vm()["bccoeff"].template as<double>();
        auto clDef = CL_DEF_MECASOL(M_time);
        if (M_weakCL)
            {
                ForEachCL( clDef,cl::dirichlet_vec,
                           form2( M_Xh, M_Xh, J )  +=
                           integrate( markedfaces(mesh,mesh->markerName(PhysicalName)), _Q<4*Order>(),
                                      //terme classique
                                      - trans((dF*Sv + Fv*dS)*N())*id(v)
                                      //terme de symetrisation
                                      - trans((dFF*Sv + Fv*dSS)*N())*idt(u)
                                      //            - trans((dF*Sv + Fv*dS)*N())*idt(u)
                                      //            - trans((F*S)*N())*idt(u)
                                      + penalisation_bc*trans(idt(u))*id(v)/hFace() ) );
            }


#if 0
        form2( M_Xh, M_Xh, J )  +=
            integrate( markedfaces(mesh,mesh->markerName("Fixe")), _Q<4*Order>(),
                       //terme classique
                       - trans((dF*Sv + Fv*dS)*N())*id(v)
                       //terme de symetrisation
                       - trans((dFF*Sv + Fv*dSS)*N())*idt(u)
                       //            - trans((dF*Sv + Fv*dS)*N())*idt(u)
                       //            - trans((F*S)*N())*idt(u)
                       + penalisation_bc*trans(idt(u))*id(v)/hFace() );
#endif
#if 0
        form2( M_Xh, M_Xh, J )  +=
            integrate( markedfaces(mesh,mesh->markerName("Force")), _Q<4*Order>(),
                       //terme classique
                       - trans((dF*Sv + Fv*dS)*N())*id(v)
                       //terme de symetrisation
                       - trans((dFF*Sv + Fv*dSS)*N())*idt(u)
                       //            - trans((dF*Sv + Fv*dS)*N())*idt(u)
                       //            - trans((F*S)*N())*idt(u)
                       + penalisation_bc*trans(idt(u))*id(v)/hFace() );
#else
        /*    form2( M_Xh, M_Xh, J )  +=
              integrate( markedfaces(mesh,mesh->markerName("Force")), _Q<4*Order>(),
              //terme classique
              - 0.0*trans((dF*Sv + Fv*dS)*N())*id(v)
              //terme de symetrisation
              - 0.0*trans((dFF*Sv + Fv*dSS)*N())*idt(u)
              //            - trans((dF*Sv + Fv*dS)*N())*idt(u)
              //            - trans((F*S)*N())*idt(u)
              + 0.0*penalisation_bc*trans(idt(u))*id(v)/hFace() );

        */
#endif

#if 1
        double Newmark_gamma=0.5;
        double Newmark_beta=1/4.;

        auto cst_acc_displ = 1./(Newmark_beta*std::pow(M_bdf_displ_struct->timeStep(),2));

        form2( M_Xh, M_Xh, J ) +=
            integrate( elements(mesh), _Q<2*Order>(),
                       cst_acc_displ*trans(idt(u))*id(v) );
#endif
        J->close();

        //  J->addMatrix(1.0, this->OpLin );

        if (!M_weakCL)
            {
                auto dirElem = this->getDirichletElement();
                ForEachCL( clDef,cl::dirichlet_vec,
                           form2( M_Xh, M_Xh, J ) +=
                           on( markedfaces(M_Xh->mesh(), M_Xh->mesh()->markerName(PhysicalName)) , *dirElem, R, vf::vec(vf::cst(0.),vf::cst(0.)) ) );

            }

        if (M_LookComment) std::cout << "[StVenantKirchhoff] : updateJacobian finsh\n";
#endif
    }

} // Feel


// instantiation
template class Feel::StVenantKirchhoff<2,2,Feel::Simplex>;
//template class Feel::StVenantKirchhoff<3,2,Feel::Simplex>;

